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ABSTEACT 


Finite Element Method has been applied to the 
Inverse and Direct Problems of Sabsonic steady and com- 
pressible now of inviscid fl\iLd in turbomachines. A 
recursive scheme for the solution of three versions of 
Indirect Problem (which have varying degree of non-linea- 
rity) and the Direct Problem has been suggested and the 
related formulations have been obtained. The Inverse 
Problem is solved on computer using a computer programme 
which has the capacity to solve all the three versions 
of the problem and good agreement of results with reported 
data is found. On the basis of the severe non-linearity 
of the governing equations for the Inverse Problem and 
the suitability of the solution scheme for the same, 
it is hoped that the Direct Problem (foraiulated in this 
work) and other similar non-linear Elliptic problems will 
also find the scheme suitable. 
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INTROLJCTION 

1,1 The Turboraachine Plow Problem 

The usefulness of obtaining solutions to turbo- 
machine flows is beyond dispute. The two main problems of 
interest are 

i) the Inverse Problem vriiich seeks the blade geometry 
for chosen velocity distribution along the blade 
stream lines and 

ii) the Direct Problem which seeks the solution to tbe 
flow through given blading and casing geometry. 

The solution of real three-dimensional unsteady 
flow has not yet been obtained even numerically let alone 
analytically. Approximate solutions have been obtained by 
patching together solutions in the meridional or hub-shroud 
plane and in a cylinderical blade to blade surface. The 
present work is concerned with the numerical solution of 
the direct and inverse problems for the blade to blade flow. 

Until recently, finite difference methods were 

exclusively used to obtain turbomachine flow solutions. 

* 

Adequate ref erences , are available in the literature for 
such methods. The ability of finite element method to deal 
with difficult boundary conditions vjith relative ease has 
attracted workers to apply such methods to the turbo- ■ 
machinery flows. References to such works are given in 
Section 1,2. 
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The present work attempts to deal with the 
Inverse and Direct problems in a unified manner* A success- 
ful computer programme has been written for the Inverse 
problem. The programme uses linear triangular elements, 
generates it's finite element mesh, uses semiband storage 
mode for storing useful stiffness matrix elements and 
solves the system of linear simultaneous equations using 
a version of Gauss elimination algorithm developed for 
sjTmnetric coefficient matrices stored in semiband mode* On 
the basis of the mathematical similarity of the equations 
governing the Direct and the Inverse problem, it has been 
anticipated that the formulation proposed for the Direct 
problem may give good results. 
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1 .2 The Finite Element Method and Compressible Flow 

The governing equations for compressible flovr 
of fluids are elliplic or hyperbolic depending on whether 
the flow is c ^ ■ .Lv.c.'ic. The equa- 

tions have varying degree of nonlinearity; therefore iiiere 
is no general method available for the solution df compre- 
ssible flow problems using the finite element method. 

Work has been reported in the field of subsonic and super- 
sonic compressible flow by many authors and the general 
direct problem of fluid flow has been solved using a 
variety of finite element formulations. As it has always 
been, more ’'rork has been done in the field of subsonic flow 
than in that of transonic flow because of additional compli- 
cations involved in the treatment of the latter. Very 
little work has been ropoibed in the direction of solution 
of the inverse probleLi of compressible fluid flow. The 
present work deals mainly witi the solution of subsonic 
cori5>ressible floi^ in a turboraachlne and therefore, in tliis 
section, only the subsonic 2-D flow has been discussed, 

M.J*. 0' Carrol and L,A. Morgan 3 have investi- 
gated iiiG utility?- of linear finite elements (Triangular 
shape) for the analysis of flow in twisted passages and 
found that faster convergence is achieved compared to 
finite difference solution for the same accuracy. 

K. Washizu and M, Ikegawa[323 have solved the integral 
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equations governing lifting surfaces in flo\^ field success- 
fully'' aiid it has hecn claimed that their method can he 
extended to other integral equations of interest. H.B. Awhi 
and J.H, owannell h.ave investigated the flow round a 
circular cylinder in a wind tunnel and found good agreement 
between the finite element and the experimental results, 

G,F, Carey I 3 Q has used a perturbation expansion to reduce 
nonlinear flo\^ equations to a sequence of linear problems 
with positive definite quadratic functionals and nonlinear 
lower-order nonhomogeneous expressions, 

T.J, Chung, J.T. Oden and S,T. 
investigated the usefulness of finite element method 
developed by Aguirre - Bamirez, Oden and V/u[^6^in solving 
the problem of unsteady flow of a rarefied gas through a 
constant area duct of irregular cross-section, S.T. Shen[^7^ 
has explored the compressible flow phenomenon from the 
view point of incorporating all the numerical requirements 
through the finite element method for the solution of the 
problems of interest. B. Davis and J.A, Hendry 1^8]]] have 
suggested the use of a global variational method with trial 
functions not satisfying the prescribed boundary conditions 
for the problem of a two-dimensional channel design ahd 
their method can find a ready finite element iirplementation. 



5 


1.3 The Inverse Problom 

The inverse problem in the context of tarbomacliines 
involves the design of blade profile for desired performance. 
The two most favoured approaches for its solution entail 
transfomiabion of the governing equations to the potential 
(w = jZ5 + i W ) plane and to the hodograph plane. 

Significant work has been reported in the former 
class by Costello Stanitz Payne 

this approach, a suitable velocity profile is chosen to 
obtain the corresponding blade profile. Costello uses the 
velocity distribution to select a suitable incompressible 
potential flow around a unit circle and then uses Lin's 
02^ transformation to transform the incompressible flow 
into linearized coii 5 )ressible flow about a cascade. Stanitz 
has used relaxation method and Green' s function to solve 
a combined continuity-irrotationality equation. The use 
of Green's function in this case has been facilitated by 
the linearization of the governing equation using tangent 
gas approximations and the results were computed at the 
channel walls, the entrance and the exit. The method was 
improved upon by Payne who u.sed a slightly different 
dependent variable without changing the nature of Stanitz' s 
linearized equation. 

The advantage of using the hodograph plane 
approach lies in the linear equations which govern the 
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potential and the stream functions in the hodograph plane. 
Significant work in this direction has heen reported hy 
Cantrell and Fowler ^3*} j Uenishi , Cohen (j^l ? an<i 

Payne|^1'l'3 • Uenislii has reported reasonable agreement 
between theor;^ and experiment for a wide variety of Dlades. 
Cohen has used constant surface velocity aerofoils as test 
cases to sliow the validity of his method; however, Payne has 
cautioned against a hasty interpretation of Cohen’s results. 



CIIAPTE?. 2 


FINITE ELEIK'IT FORt''IULATION OF INVERSE PROBLEM 


Tlis govGrriing equcition for tliG ■bwo—d.iinGnsional 
comp possible, irrotational, inviscid and steady flow 
through, a channel can be written as follows 5 [-101 = 

. 22 .. . 2 . 2 


where 


3 


(2.1) 


A = (1 + ( / (1 - (- 


, -(k + 1 / k - f 

025 


B = . c2 (1 + C^4-l) c2) / (1 - (^-= 5 !) o2) 


(2 k /k - 1) 


C = - / (1 - ( 


k " 1 


■) C'^) 


The domain of the problem is a rectangle in liie 

el - plane which has been subdivided into three-noded 
^ ^ S' In C 

triangular elements. The values of A, B, C, 

^In- C arc evaluated at the nodal points using the aval- 

lable taiovm values of C and thus the governing equation is 

rendered a quasl-polsson equation. To evaluate tte quan- 

1 B c ^ ^ and ^ - ■ at the nodal points, 

titles A, B, o , ^ ^ ^ 

the following shape function is chosen to approximate the 
behaviour of C in an element J ^ 

In G iy Ug = 3C ^il 


( 2 . 2 ) 



8 


where 

K1 =[^^32 < 0 - 02> - % - V'aO /2^e 

N2 = 0^31 ( 0 - 03) + 031 ( - 'P3) 3 / 2 

N3 =0^21 ^ ^ " ^1^ “ ^21 / 2 A^ 

and Uj^ are nodal values of U, 


The variational fimctional corresponding to 
Poisson equation of this form is given by 

( 


» ■ lie Hv< 




"ST 






+ ^ne °e3 <1 0 (2.3) 

where A is the value of A at the centroid of element 
— ec. — 

and D is the value of ri^t hand side terms of equation 
ne 

(2,1) at the nodes 1, 2 and 3 of the element. 


Minimization of the fimctional I gives us the 
following stiffness equation: 


C^elt 



(2. if) 


M = ^eo 1 

i 0^2 ^32 ^31 ^32 ^2l] 


^32 ^32,^1 f32.f3^ 


^31 %1 ^21 

+ 

'Psi '^31.^21 


\ Sym 02^ 


Sym. 2 ^ 


^ . J 

i 

— J 


(2.ifa) 
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0 . 


-ne 



(032^2 "^32 ^2^ %2jj 0 ^0 032{jV^‘^0 

2 ^^31 % " 031*^3^ •+• ( 03 ijjs^d 0 df-c^3^jj 0 d0 dy^)/A 

3 ^^21^1 “^21 ^1^ *^^21]] ^ 02lJJ'f''^^ dt^VA 


Di = 


^i{ “I + (li + 1)/^j_ / 2^<<^ /~l^1 - (k - '>)/6i/2j- 
-^i /{ ^ / 2]<2 ? 

^ ^ y^32 ^1 ”^31 ^2 '*'^21 ^3 ^ ’ 


^2 

^1 

2 


Ji 

'ec 

U_ 


0 


^ ^32 ^1 “ ^31 ^2 ^21 ^3 ^ ’ 

^1 ^^32 ^ 21 *" 1^32 ^21 ^ ^ ’ 

e 

^2 ^^ 32 !^ 31 “ ^32 ^31 ^ \ 5 

U3 (031 IP 21 "ip31 021^ / ■^e ’ 

e 

r 


(2.4Td) 
2k/k-1 

(2.5) 

(2.5a) 

(2,5b) 


d0 df = 0c • ^e’ 


J'pc30 d4/ = . A^ 


2 TJ. 


= (1 + 0.5 (k + 1) e ^) / (1 - 0.5 (k - 1) e 


(2.5d) 
(2.5e) 

(2.6) 

2U (k+l)/(k-1 
"" ) 


= U 


0 = 00 i V = V' ^ ’ 


and 0^ = (0^ +02+03) / 3 ; 1^3 = (V^l + V'a +^^ 3 ) / 3 


C 2 . 7 ) 



CHAPTER 3 


solution technique 

Some variations of the forniiiLation presented 
in the previous chapter are used to solve the following 
three versions of the governing eauations s 


1 ) 


>2 


In C 


9 2 


00*^ 




In C 

.2 


0 


(3.1) 


2 ) 


where k = - ^ 


and 


.5^ 


1 + / 


1 + C‘ 


(Tangent gas approximation and transformation used) 


"5 ^ In C 3 ^ In 0 


^ j6‘ 


2 Vf' 


1 + c' 


f , 3 In C ^ , 3 In C s 


2>[l/ 
(3.2) 


(Only k = - 1 nsed) 


3) A P 

- 3 0^ 




3 Ip 


~WT 


3>ln G 
(3.3) 


llie following steps are talien to solve iiie above 


equations : 
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1 . Divison of the domain into triangular elements uith 

aspect ratio nearly equal to one. This is done in the 
present \rork by computer (using an automated mesh genera- 
tion scheme). 


2. Calculo-tion of clementvrisc node numbers, nodal coordi- 
nates and value of A at the element centroid in case of 
equation (3.3). 

3. Calculation of element stiffness matrix using equation 
(^-,1). The value of is taken as unity for the first 
two cases. 


4-. Assemblage of the overall stiffness matrix. 

5 . Calculation of the value of light hand side terms in 
the case of equations (3.2) and (3.3) and using this as 
Di to evaluate the load vector. The load vector is null 
in the case of equation (3.1). 


6, Assemblage of the overall load vector, 

7. Application of boundary conditions in a suitable fashion. 
The technique followed in the present worl: for incorpora- 
ting the boundary conditions is as f ollov^s i 

i) Multiply the diagonal terra of the overall stiff- 
ness matrix corresponding to the degree of free- 
dom which has to be restrained by a very large 
10 

number, say 10 . 

ii) M-ultiply the prescribed boundary condition with 


the corresponding modified diagonal blfement '^" 









tlie overall stiffness! matrix and replace the 
respective element of the load vector hy tiiis 
number. 
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8, Solution of the obtained sot of linear simultaneous 

equations using a suitable scheme. In tiie present work 
a Gauss elimination equation solving routine has been 
developed for symmetric banded coefficient matrices 
stored in semiband form, 

9* Comparison of the currently obtained solution with the 
previous solution to find if the solution has converged. 
This is accomplished by compa.r:ng the cumulative relative 
difference of nodal values with chosen permissible error. 

10. Stopping the iterative process if convergance is achieved. 
Otheri'fise, rename the new solution as old solution and 
repeat steps tvo through ten. 



As mentioned above, the cumulative relative 
difference is con^ared with a chosen permissible error para^ 
meter. The cumulative relative difference is calculated 
as follows ; 


GUP. (Cumulative relative in^rovement parameter) 


In “ In C^T , 

new old 

In C 


where N is the total number of nodes 



CHAPTER h 


RESULTS Al® CONCLUSIONS 

The results obtained for versions 1, 2 and 3 are 
presented in Tables 1, 2 and 3 respectively. Table 4 shows 
the values of C,R.I,P. and the central processor unit time 
required per iteration for versions 2 and 3* Figure t 
depicts the solutions of all the three versions. 

It is seen that the method used gives almost 
identical results for versions 1 (Eqn. 3*1) and 2 (Eqn. 3*2) 
of tile inverse problem. It is also seen that the use of the 
same interpolating function to solve the third version of 

the problem gives results which compare very well with the 

♦ 

results reported by Stanitz|^10^ . Ihe solutions of versions 
2 (Eqn, 3«2) and 3 (Eqn. 3*3) are found to converge after 
two itei*ations only when the convergence criterion is average 
one percent improvement in the nodal values and the order of 
magnitude of C.R.I.P becomes unity on the con^jletion of the 
second iteration. 

It is found that no improvement of values is 

effected by carrying out more than six and seven iterations 

in the case of versions 2 and 3 respectively and in both 

*.3 

cases the order of magnitude of C,E,I,P, becomes 10 on the 
completion of the third iteration. 


* ttOiSam Mffme at the aest 
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After the stud^r of tlie resifLts obtained, the 

following is concluded; 

1 . The same degree of accurac 3 ?' can be obtained by solving 
versions 1 and 2 using the same typo of elements and 
j.nterpolating function. This renders the analytical 
manipulations based on the tangent gas approximation 

as done by Stanitz Qo] for the linearization of the Eqns. 
(i.e,, getting version 1) unnecessary, particularly in 
view of the very siiiall number of iterations required for 
the convergence. 

2. The version 3 is a better choice in comparison v;ith the 
versions 1 and 2, The excess con^jutational time needed 
in the case of the former is not much and also this 
excess tirae is more than compensated by the ingjrovement 
obtained in the nodal values. 

3. Further in^^rovement in the values obtained can be 
effected by using a higher order element and an itera- 
tive scheme for solution of linear equations . The 
latter may require special iterative techniques suited 
to the semi-band storage of the symmetric and banded 
stiffness matrix. Another manner in which the in 5 )rove- 
ment can be expected is through the use of quadrilateral 
element in place of the linear triangular element. 

As far as the present problem is concerned the conver- 
gence requirement of one percent improvement per node 
(on the average) seems to be quite adequate. 
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TABLE 4 


lATIC 

fo. 

LY E R s i 0. 

N TWO 

0 VERSION 

THREE 


1 C.P.U 

] C.B.I.P 1 

C.P .U 
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ap?si:dix a 


A.1 PiiO Direct Problo?:: 


n-Tic velocity potentlaJ- function 0 satisfies the 
1 'ollowir.y partial dili'crential equation for tv^o- dimen s ional , 
steady, inviscid compressible flov/s : 


0 •¥ 0 
^x>: ^yy 


^yy ^ ^ ^y ^xyj 


where 


= C5^- 


2 2 

1 (k - 1) (0^^ + 0y) 


(A.1.1) 


Tlierc is no universa.lly valid technique available 
as yet to solve the above equation. Gostellow has reviewed 
the various classes of this problem. In this reference, 
efforts have been directed towards the performance analysis 
of medium solidity cascades termed as the general direct 
problem. There are four main categories of solution; namely, 
series solution methods, iterative methods, matrix methods 
and streamline curvature methods. 

Two classes of work utilizing series solution 
methods have been reported out of which one treats the velo- 
city potential function as the dependent variable while the 
other uses the streamline function for the same purpose. 

The dependoit variable is approximated by a series of the 
following type 
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+ 1 
n 


-2n 


(A.1.2) 


4 m 


^ « • « 


C. 


• ■ if; eh-.’.ont variable jZ5 or ^ and M is 
r . Icr, is the incompressible solution 
.',ru i’-Lncticns of the lover order coefficients 
are the independent variables. 


'l’:c si; ul''.tion of the right hand side of the 
iavci’::in, ir'nntlcn V^' a source - sink distribution within 
tb.c itel, of inco!.r,5rossiblc flow around the profile and the 
calc .i.-iLion t’'x*rcof in iterative steps is the essence of the 
it.orative * cthods vniich arc reported to converge slowly as 
the non! ' corulitior* is approached. 

In ijost matrir. methods, the governing equation is 
solved ar, if it vere a poisson equation using ten~point-star 
finite difference scheme or relaxation. The various formu- 
lations differ in their treatment of density gradients. 

The streamline curvature methods use streamlines 
and ’Quasi-orthogonals’ (lines passing from one channel wall 
to the otlier) . The input data is differentiated twice and 
the velocities obtained are then used in the velocity gra- 
dient equation. A start is made v;ith approximate velocities, 
along the midstream line obtained from the previous iteration . 
A numerical integration is then performed in both directions 
to find the blade surface velocities and then the velocity 
levels are adjusted to satisfy tlie continirity equation# 
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tc 


• Kt. . 1,,^ .c methods have also been applied 

■; ; i.eir Vriof description appears in 


• cc. 


• .> * V 


il u': 'c:. '^n;: Soj-ition Technique 


vvcridn/ equation is rewritten as follows; 




v/h.orc 


^xx ^yy ^ ^y ^xy 1 


and 


Co‘ 


- 2 2 

Tj (h •« 1) (0„ + 0„) 

c. JL Y 


(A. 2.1) 

The flow field around tiie profile is divided into 
straiclit sided triangular elements and the behaviour of the 
function 0 is approximated with a quintic polynomial in 
the aroa coordinates Z^v ^2 ^3’ 




a. 2.2) 


where are coefficients of the terms of the quintic 

polynomial which are evaluated in terms of the nodal values 
of the degrees of freedom as follows: 

The values of 0 as referred to the physical 
X « y coordinate system can be evaluated from „ 7 

2 

in the following manner: 


(A.2.3) 
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0e 


= ("if 


(A.2.^) 

(A.2.5) 
(A. 2, 6) 


As before, the value of L for an element is 
confuted from the known nodal values of 0 and the usual 
variational-finite element technique is used to get the 
stiffness equation as follows: 


where^ 


0 

0 


,Z2 

= i 

Qi\ 

. Z^, 

^2 


(A. 2 . 7 ) 
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^0^ 

3^ 02 

3^01 
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5^1 ^2 ^ 

9 z| 


3 02 

3 02 

3^02 

3^ 02 

0 

2» 

3 2i 


’ 9 ’ 

3^1 ^2 

’ '^4 


S 03 

a 03 

9 ^ 0:^ 

3 ^ 03 

II2 

‘3» 


■ ' 3 Z 3 

’ 3 zf 

92 ^ Z 2 ’ 

9Z| " 





T 

(A. 2.8) 

'lk ■ 

il » 

L^2k*l» 

•L.' 03k i 

. ! 




k = 

= 1. 2. 3, 

if & 5 

(A.2,9) 


K] 


here subscripts 1, 2, 3 and if refer to the three 
nodes of the straight sided triangle and its centroid 
re sp ectiv ely • 
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^ = - F (Z^ , Zg) dA 

K 

ff-, ® 2A.nl mill 

I Z^ Z| Z” dA = ' (1% m 1- n H- 2l - 


(a.2.10) 


(A.2.11) 


-3 f 


^ r Sf "Sf 


TL 


ax 2)y 


IdAlCs'l] 




(A. 2. 12) 


where 


P]- ifCft «3' D 


^ f Bf 

,ax -5 7 


1 


■ ff c,M ^1’ r'-i’r’x*?? 


(A.2.13) 


where CT is the well known Jacohian matrix snd. 


0. . ^ 

0 h J 7 7 

j ‘j-j }^2 


^32 ^13 


■ L' l?i il 


hifere j =1,2, 3s ^ -1}2, 35^>? 


(A.2.14.) 


(A.2.15) 


T = 


2 ^13 ^13 


x.o Xo7f^l3 ^23 


‘13 23 


?13 ^23 
^ ^23 ^23 


A. 2. 16) 


^13 ^23 
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The 0^ term suppressed soon after computation of element 
stiffness matrix. 


f (Z^ , 22? ^3^ “ f (Z^ , Z^) 


=: 1 > , 

^2» ^2’ 

Zg Z3, z 

3» ^2 - ^2 

^2 ^3 

- Z3 Zg, 

Z3 - zfz,, 

Z^ Z^ 2^, z^ Z2 

- 2^ 

z| Z3 - z^ 

z 3 , z 3 Z3 

7 73 72 ™2 

- z.^, ^2 ? 

7 ^ Z^ 
Zp ^3 

72 1-2 y 3 

, z^, 

z| - z2 z|, 

„3 72 „2 „2 

^2 ^3 ■“ ^^2 ^3 ’ 


- Z^ z| (A, 2. 18) 



The quasi-Poisson equation (A,24l) is solved 
iteratively using a suitable method for the solution of -the 
generated set of linear simultaneous equations in the follow- 


ing manner? 

C‘]{'r 



CA.3.1) 



3 ^ 


Stops taken to effect this are as follows J 

1. Computation of element stiffness matrices. 

2. Assembly of overall stiffness matrix and its storage in 
the disk, 

3. Calculation of element load vector on the basis of 
the current values of 0^^ & its assembly. 

Retrieval of overall stiffness matrix from disk for ite- 
ration numbers greater than one and/or application of 
the boundary conditiCns besides enforcing the_e<5uality of 
boundary conditions v/herever prescribed. 

5. Solution of the system of linear simultaneous equations. 

6. Testing of the new solution using a suitable convergence 
criterion. 

7. stopping the process if the convergence criterion are 
fulfilled otherwise repetition of steps 3 through 7. 

Convergence Criterion 

The convergence parameter is compared with the 

cumulative relative improvement of nodal values in the same 

way as has been done in the case of the inverse problem. 
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APPENDIX B 


PBOGRAMME LISTING AND USERS’ 
INSTHJGTIONS 


B.1 TT .c; AT R « Inst ructi ons. 


This ppogrs-inin© has hssn Usvslopsd. for rinning 
on DEC 10 Time sharing system and it \ises the disk storage 
for retrieving and filing data. The various inputs to the 
programme have the following format and significances 

^ ^ jjp integer number. It controls the flow of 

computations in respect of different versions 
of the non-linear problem. 


UP = 1 Linearised Problem (Version 1 ) 

UP 2 Partially Linearised Problem (Version 2) 

UP = 3 Complete Non-linear Problem (Version 3) 

UIT Do loop parameter which fixes the upper limit 
on the n-umber of iterations. 


NH 

NV 


PHII & 
PSIF 


Number of nodes on the horizontal side of the 
domain. 

Number of nodes on the vertical side of the 
domain. NH is greater than NV if not equal 
to the latter. 

Range of values of independent variable 
represented on the vertical side. 
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DEL 


Q(I), 1=1,11 


XP 

The data is 

aUBROUTINES 
MAIN 1 and 
MAIN 2 
MAIN 3 
PLACE 


Danping factor* Taken Mnitj if damping is 
not needed. 

Tobal nui’abcr of Boundary nodes. Calculated 
by the Computer itself • 

Vector of length N. Having the values of 
the prescribed Boundary conditions as it's 
elements. Values are read in the following 
way: 

1, Nodal values at the nodes on the left 
hand side 

2, Nodal values at the nodes on the right 
hand side 

3, Top side nodes leaving those already 
accounted for from left to right. 

if. Bottom side nodes in the sarae manner as 
in case of 3 above. 

The limiting permissible cumulative relative 
error (In subroutine CONVER) 

read from a file named DAT^f CDR from disk. 

Compute the nodal numbers and 
coordinates . 

Compute Ag^ for NP = 3 (i.e., Version 3 ). 
Places the boundary conditions appropriately 
from vector Q read in the main programme 
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SOLVE 


CONVER 


WRITE 


into the vector PHI and generates a vector 
NO of length N whose elements indicate the 
nnmbor of the nodes on the boundary. 

Uses a modified Gauss Elimination algorithm 
suitable for syrametric matrices stored in 
semi-band mode to solve tie stiffness equa- 
tion . 

Tests convergence and renames the new solu- 
tion as the old solution before returning 
to the main programme. 

Writes the results in a file named DESIGN.DAT 
and stores in users* area. 
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GOTO 6 
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RETURN 

END 
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CUMM0N/B1/PHU434) ,0r434),nR(434,9),NUCl3S),NN0{3) 
rt3MMO,^/tt2/Ntt,NV,N^,i^F,NVE,TLT,N,THr>i«l,irot.,DX,DY»XC,TC,BKC3,37,Xt33 
UY(33, 01(3), B(35, DEI 
s> TxPE 98 

98 FORMA'H' enters CULVER') 
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